About this code

#load libraries
library(raster)
## Loading required package: sp
library(rhdf5)
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

Let’s Open a File!!

f<-"../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"
# view h5 structure
h5ls(f)
##    group                                 name       otype  dclass
## 0      /                     ATCOR_Input_File H5I_DATASET  STRING
## 1      /                 ATCOR_Processing_Log H5I_DATASET  STRING
## 2      /                Aerosol Optical Depth H5I_DATASET INTEGER
## 3      /                               Aspect H5I_DATASET   FLOAT
## 4      /                          Cast Shadow H5I_DATASET INTEGER
## 5      / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6      /                 Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7      /                  Illumination Factor H5I_DATASET INTEGER
## 8      /                          Path Length H5I_DATASET   FLOAT
## 9      /                   Processing Version H5I_DATASET  STRING
## 10     /                          Reflectance H5I_DATASET INTEGER
## 11     /                Shadow_Processing_Log H5I_DATASET  STRING
## 12     /                      Sky View Factor H5I_DATASET INTEGER
## 13     /               Skyview_Processing_Log H5I_DATASET  STRING
## 14     /                                Slope H5I_DATASET   FLOAT
## 15     /          Slope_Aspect_Processing_Log H5I_DATASET  STRING
## 16     /                  Solar Azimuth Angle H5I_DATASET   FLOAT
## 17     /                   Solar Zenith Angle H5I_DATASET   FLOAT
## 18     /                    Surface Elevation H5I_DATASET   FLOAT
## 19     /                 Visibility Index Map H5I_DATASET INTEGER
## 20     /                   Water Vapor Column H5I_DATASET INTEGER
## 21     /             coordinate system string H5I_DATASET  STRING
## 22     /                       flightAltitude H5I_DATASET   FLOAT
## 23     /                        flightHeading H5I_DATASET   FLOAT
## 24     /                           flightTime H5I_DATASET   FLOAT
## 25     /                                 fwhm H5I_DATASET   FLOAT
## 26     /                             map info H5I_DATASET  STRING
## 27     /              to-sensor azimuth angle H5I_DATASET   FLOAT
## 28     /               to-sensor zenith angle H5I_DATASET   FLOAT
## 29     /                           wavelength H5I_DATASET   FLOAT
##                dim
## 0                1
## 1                1
## 2    544 x 578 x 1
## 3    544 x 578 x 1
## 4    544 x 578 x 1
## 5    544 x 578 x 1
## 6    544 x 578 x 1
## 7    544 x 578 x 1
## 8    544 x 578 x 1
## 9                1
## 10 544 x 578 x 426
## 11               1
## 12   544 x 578 x 1
## 13               1
## 14   544 x 578 x 1
## 15               1
## 16           1 x 1
## 17           1 x 1
## 18   544 x 578 x 1
## 19   544 x 578 x 1
## 20   544 x 578 x 1
## 21               1
## 22         5732053
## 23         5732053
## 24         5732053
## 25         426 x 1
## 26               1
## 27   544 x 578 x 1
## 28   544 x 578 x 1
## 29         426 x 1

Import Spatial Information

# import spatial info
mapInfo <- h5read(f,
                  "map info", 
                  read.attributes = TRUE)
mapInfo
## [1] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## attr(,"Description")
## [1] "Basic Map information for envi style programs"

Grab Reflectance Materials

## read in reflectance data attributes

reflInfo <- h5readAttributes(file = f, 
                             name = "Reflectance")
## define scale factor
scaleFactor <- reflInfo$`Scale Factor`

## define no data value
noDataValue <- as.numeric(reflInfo$`data ignore value`)
str(noDataValue)
##  num 15000

Import Data Dims

# Open the file for viewing
fid <- H5Fopen(f)
# open the reflectance data
did <- H5Dopen(fid, "Reflectance")
# grab the dataset dimensions (column, row, band)
sid <- H5Dget_space(did)
dims <- H5Sget_simple_extent_dims(sid)$size
dims
## [1] 544 578 426
# close all open connections
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)

Read in Reflectance Data

# extract slce of H5 file
b56 <-  h5read(f,
               "Reflectance",
               index = list(1:dims[1],1:dims[2],56))

Convert to Matrix

# convert to matrix
b56 <- b56[,,1]

Plot Data

## let's plot some data

## first change from scientific notation to decimals
options("scipen"=100, "digits"=4)

## now try an image
image(b56)

## log transform it--too dark, can't see well otherwise
image(log(b56),main="log transformed data")

## histogram
hist(b56,
     col="springgreen",
     main="Distribution of Reflectance Values \nBand 56")

## Data Clean Up

# assign no data value to object (already created noDataValue earlier)
# set all reflectance values = 15,000 to NA
b56[b56 == noDataValue] <- NA

# apply scale factor (created earlier)
# reflectance values between 0-1 (the valid range)
b56 <- b56/scaleFactor

# view distribution of reflectance values
hist(b56,
     main="Distribution with NoData Value Considered\nData Scaled")

Transposing (flipping) the data

# Because the data import as column, row but we require row, column in R, we need to transpose x and y values in order for our final image to plot properly
b56 <- t(b56)
image (log (b56), 
      main="Band 56\nTransposed Values")

## Define spatial extent

# We can extract the upper left-hand corner coordinates.
# the numbers as position 4 and 5 are the UPPER LEFT CORNER (x,y)
mapInfo <- strsplit(mapInfo, ",")
mapInfo <- unlist(mapInfo)
mapInfo
##  [1] "UTM"               "1"                 "1"                
##  [4] "325963.0"          "4103482.0"         "1.0000000000e+000"
##  [7] "1.0000000000e+000" "11"                "North"            
## [10] "WGS-84"            "units=Meters"
# grab the X,Y left corner coordinate ensure the format is numeric with as.numeric()
xMin <- as.numeric(mapInfo[4])
yMax <- as.numeric(mapInfo[5])

# we can get the x and y resolution from this string too
res <- c(mapInfo[2],mapInfo[3])
res <- as.numeric(res)

# finally calculate the xMax value and the yMin value from the dimensions
# we grabbed above. The xMax is the left corner + number of columns* resolution.
xMax <- xMin + (dims[1]*res[1])
yMin <- yMax - (dims[2]*res[2])

# also note that x and y res are the same (1 meter) Now, define the raster extent define the extent (left, right, top, bottom)
rasExt <- extent(xMin, xMax,yMin,yMax)

# now we can create a raster and assign it its spatial extent
## CRS is a coordinate reference system
b56r <- raster(b56,
               crs=CRS("+init=epsg:32611"))
# assign CRS
extent(b56r) <- rasExt

# view raster object attributes
b56r
## class       : RasterLayer 
## dimensions  : 578, 544, 314432  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326507, 4102904, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0033, 0.5678  (min, max)
# plot the new image
plot(b56r, main="Raster for Lower Teakettle \nBand 56")

Import NEON Functions

library(devtools)
#install_github("lwasser/neon-aop-package/neonAOP")
library(neonAOP)
??read_band
??open_band

b55 <- open_band(f, 
                 band = 55,
                 epsg = 32611)
b55
## class       : RasterLayer 
## dimensions  : 578, 544, 314432  (nrow, ncol, ncell)
## resolution  : 1, 1  (x, y)
## extent      : 325963, 326507, 4102904, 4103482  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0029, 0.5622  (min, max)
## plot data
plot(b55, main="Raster for Lower Teakettle Band 55")

# import several bands
bands <- c(58, 34, 19)

## create a raster stack
RGBstack <- create_stack(f, 
                         bands = bands, 
                         epsg = 32611)
plot(RGBstack)

# plot RGB stack
plotRGB(RGBstack,
        stretch="lin")

#cir image
bandsCIR <- c(90, 34, 19)
## create a raster stack
CIRstack <- create_stack(f, 
                         bands = bandsCIR, 
                         epsg = 32611)
plot(CIRstack)

# plot RGB stack
plotRGB(CIRstack,
        stretch="lin")